function [eta, gamma, bMarkup, cost] = rcnl_eta( ...
    kappa,   ...
    cdindex, ...
    cdid,    ...
    firmid,  ...
    prodid,  ...
    yearid,  ...
    bmc,     ...
    s_jt,    ...
    p_jt,    ...
    w,       ...
    derMat,  ...
    fes      ...
)

bc = (yearid(cdindex) > 4) .* kappa;
bMarkup = zeros(size(s_jt));
for m = 1:max(cdid)
    sel = (cdid == m);

    ownmat = repmat(firmid(sel), 1, sum(sel));
    owner  = ownmat == ownmat';
    omega  = owner + (1 - owner) .* (bc(m) * (bmc(sel)) * (bmc(sel))');
    der    = derMat(prodid(sel), prodid(sel), m);

    bMarkup(sel) = - (omega .* der') \ s_jt(sel);
end
cost = p_jt - bMarkup;

% Linear parameters
nfe    = size(fes, 2);
xx     = fes' * fes;
wtilda = w(:, 1:end-nfe) - fes* (xx \ (fes' * w(:, 1:end-nfe)));
ctilda = cost - fes * (xx \ (fes' * cost));
gamma1 = inv(wtilda' * wtilda) * wtilda' * ctilda;

% Fixed effects
y      = cost - w(:, 1:end-nfe) * gamma1;
gamma2 = inv(xx) * fes' * y;
gamma  = [gamma1; gamma2];

% Unobserved cost
eta = cost - w * gamma;
